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Abstract 


This report describes the preliminary results of an investigation on component reliability analysis 
and reliability-based design optimization of thin-walled circular composite cylinders with average 
diameter and average length of 15 inches. Structural reliability is based on axial buckling strength 
of the cylinder. Both Monte Carlo simulation and First Order Reliability Method are considered 
for reliability analysis with the latter incorporated into the reliability-based structural optimization 
problem. To improve the efficiency of reliability sensitivity analysis and design optimization 
solution, the buckling strength of the cylinder is estimated using a second-order response surface 
model. The sensitivity of the reliability index with respect to the mean and standard deviation of 
each random variable is calc dated and compared. The reliability index is found to be extremely 
sensitive to the applied load and elastic modulus of the material in the fiber direction. The cylinder 
diameter was found to have the third highest impact on the reliability index. Also the uncertainty in 
the applied load, captured b\ examining different values for its coefficient of variation, is found to 
have a large influence on cylinder reliability. The optimization problem for minimum weight is 
solved subject to a design constraint on element reliability index. The methodology, solution 
procedure and optimization esults are included in this report. 

Shell Buckling Analysis 

The circular composite cylir ders considered in this study tall under the general category of thin 
shell structures with the displacement field described by the first-order shear deformation theory 
formulated as 

u(x,y,z ) = MoU.y) + z< Px( x - v) 

v(x,y,z) = v 0 (x,y) + z<py(x,y) ^ 

w(jc,y,z) = w 0 (x,y) 


where uq, vq, wq represent the midplane displacements in x, y, z directions (see Fig. 1), respectively, 
and <j> x ,<(>y describe rotations about the y and x axes, respectively. 

-N x 


Edge 3 
-Nx 


Fig. 1 Computational model used for cylinder buckling analysis 

The strain-displacement relations are based on Sanders-Koiter shell theory with the in-plane strains 
(£ x ,£y,Y%’ at z- 0), transverse shear strains (Yxz'Yyz' dlz = 0), and curvatures (K x ,K y ,K xy ) 
formulated as 
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where R is the shell radius of curvature in y direction (see Fig. 1). Hence, the elastic strain energy 
stored in the shell under axial compression can be expressed as 
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where Ay,By,Dy are the laminate extensional, coupling, and bending stiffness matrices, 
respectively, and C pq is the transverse shear stiffness matrix with the strain vector defined as 


{e} 7 = {e? 4 Y°xy K x K y K xy rfz xjz} 

The work done by the applied edge load shown in Fig. 1 can be expressed as 



(4) 


(5) 


where N x is the stress resultant in the x direction. 

To obtain the critical buckling load, the displacements mo,vq, vv 0 and rotations <p x , <p y are 

approximated by different Rilz series with Legendre polynomials used as the interpolation 

functions such that the essential boundary conditions are satisfied. Then by applying the principle 
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of minimum total potential energy, the critical buckling load is found by setting Eqs. (3) and (5) 
equal to each other and solving the resulting eigenvalue problem for the critical load factor such that 

N = A N ( 6 ) 

iv cr /l cr iY x 

A computer code (hereinafte r referred to as the shell code) based on the described analysis 
procedure 1 is used to calculate the axial buckling loads for circular cylinders with clamped loaded 
edges. In this case, edge shortening is allowed along edge 4 with edge 2 kept fixed (see Fig. 1). 
Also the conditions of symmetry along edges 1 and 3 are enforced by setting the v displacement 

and <j> y rotation to zero. Tabic 1 below shows the list of boundary conditions. 

Table 1 . D escription of boundary conditions for the shell m odel a 


Paiameter Edge 1 Edge 2 Edge 3 Edge 4 


u 

^ 

0 

1 

0 

0 

V 

1 

1 

1 

1 

w 

0 

1 

0 

1 


0 

1 

0 

1 


1 

1 

1 

1 

a 0 = free, 1 = 

fixed 





For validation purposes, we compared the results found using the shell code with those reported 
previously by Waters 2 as we ll as those we obtained using the STAGS structural analysis code. 
Table 2 below shows the geometric and material properties used for each specimen while table 3 
gives the laminate ply distribution and the corresponding cylinder buckling force. 


Table 2. Gee metric and material properties of each cylinder specimen 


Specimen 

L, in 

D,, in 

V in 

E,, psi 

psi 

V 12 

G 12 , psi 

1 

14 

15.75 

0.005 

18.5111e6 

1.64e6 

0.3 

0.8706e6 

2 

14 

15.75 

0.005 

18.5780e6 

1.64e6 

0.3 

0.8737e6 

3 

14 

15.75 

0.005 

18.6705e6 

1.64e6 

0.3 

0.8780e6 

4 

14 

15.75 

0.005 

19.2588e6 

1.64e6 

0.3 

0.9057e6 

5 

14 

15.75 

0.005 

18.6154e6 

1.64e6 

0.3 

0.8754e6 


The results obtained from the shell code are in fairly good agreement with those found using 
STAGS. The difference ranges from the lowest (0.3%) for specimen 5 to the highest (17%) to 
specimen 1. The large diffe rence observed in the case of specimen 1 may be attributed to using a 
Legendre polynomial that is not of sufficient accuracy for the selected ply pattern, e a so 
observed some differences with the results found by Waters . Of particular interest is the 
difference between specimens 4 and 5. We found cylinder 5 to be stronger of the two whereas the 
results obtained by Waters shows the opposite. Upon closer examination of both Si AUS ana 
shell code models, we did rot find anything that would explain the reason for this discrepancy. 
However, it appears that for the case when the cylinder length and diameter are of nearly equal 
values, the laminate pattern reinforcing the hoop direction (specimen 5) provides a suffer design 
than that reinforcing the ax al direction (specimen 4). This response is attributed mainly to the 
Poisson’s effect. Our results are consistent with the nonlinear static and nonlinear transient 
buckling analyses conducted by Hillburger and Starnes . 
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Table 3. Comparison of predicted and measured buckling loads 


Axial Buckling Force, lb 


Specimen 

Ply Distribution 

Linear 

FEA a 

Non-Linear 

FEA a 

Experiment 3 

STAGS b 

Shell 

Code 

1 

(±45 / 0 / 90) s 

44,466 

39,670 

30,164 

42,355 

49,553 

2 

(±45 / +45) 2i 

105,829 

95,618 

73,975 

104,044 

111,349 

3 

(±45/ 0/ 90), s 

158,258 

- 

147,759 

180,443 

185,420 

4 

(±45 / 0 4 / +4‘») s 

141,332 

128,610 

125,542 

154,655 

158,319 

5 

(±45/90 4 /+45) s 

139,411 

97,047 

91,667 

167,175 

167,717 


‘Values reported by Waters^” using L = 16 in. with a 1 in. cap at each end. 
b 51 quadrilateral elements along the length andl69 along the circumference. 


Structural Reliability Analysis 

The limit-state function for the critical axial buckling force is formulated as 

g(X) =P cr -P a <7> 

where P a is the resultant axial force (Load) acting on the cylinder, P rr is the corresponding buckling 

force (Resistance), and X = \X { ,X 2 X„} 7 is the vector of continuous random variables affecting 

the buckling load. According to Eq. (7), g < 0 represents failure, g>0 indicates safety, and g = 0 
represents the limit-state surf ace separating the failure and safe regions. 

The probability of failure is expressed mathematically as Pf = P(g(X) < 0), and found as 

Pf = Jf-J Q /x(*l.*2 x,,)dx { dx 2 ...dx n ( g ) 

where / X (x) is the joint probability density function of n random variables and £2 is the failure 
region defined by g < 0. 

Since the determination of the joint probability density function is not possible in this case, several 
alternative techniques may tv used to estimate the probability of failure as defined by Eq. (8). Next 
we discuss two such techniques used in this investigation. 

Monte Carlo Simulation 

To estimate the failure probability using the Monte Carlo simulation (MCS) techniques, a large 
population of response samples must be obtained. Random sampling is done by using a random 
number generator to obtain tue values for Xj, X 2 , ...,X n based on the specified mean, variance, and 
distribution type of each random variable. Using each random set as input, the shell code is used to 
obtain the corresponding buckling load followed by the limit-state function in Eq. (7) at each 
simulation cycle. Hence, the failure probability and its corresponding coefficient of variation are 

found as 


P f = N f /N 


(9) 
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COV(Pf) = 


- Md 



( 10 ) 


where N is the total number ;>f simulation cycles, and Nj is the number of cycles resulting in 
“failure” (g<0). The estimates obtained for P f and COV(P f ) become exact as N approaches 

infinity. 

In this study, we treated the applied load, geometric dimensions of the cylinder, and engineering 
properties of the material as random with statistical properties as defined in Table 4. The material 
properties correspond to carbon/epoxy (AS4 12k/3502) unidirectional tape with nominal 
engineering properties obtained from MIL-HDBK-17-2E. Due to unavailability of data, the 
coefficients of variation in the applied load and geometric variables are assumed. All random 
variables are also assumed to be statically independent. The shell laminate contains 16 layers each 
defined by separate thickness and ply angle variables. Therefore, Table 4 describes a total of 39 
random variables. 


Table 4. Statist ic i l properties of random variables affecting cylinder reliability 


Random Distribution Mean 

Variable (No) Type 


COV 

(%) 


£,( 1 ) 
EA 2) 


12 


( 3 ) 

G n (4) 

' ply (5-20) 
C(21-36) 

D( 37) 

L(38) 

PA 39) 

’Standard Deviai ion 


Normal 

18.0e6 psi 

3.19 

Normal 

1 .35e6 psi 

4.26 

Normal 

0.226 

5 

Normal 

0.543e6 psi 

5.16 

Normal 

0.005 in. 

1 

Normal 

[±45/90 4 / +45/±45/90 4 / +45] 

r 

Normal 

15 in. 

i 

Normal 

14 in. 

i 

Normal 

14.369e4 lb 

5 


Of a total of 5,3 14 simulation cycles conducted, buckling failure (i.e., g < 0) was observed in 988 
samples. This gives a probability of failure P f =0.186 for this specimen based on the assumed 
mean and variance for the resultant compressive axial load. The plots of P f and COV are shown in 
Fig. 2 with the coefficient of variation conversing to COV {Pf) = 0.0287 . 
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Fig. 2 MCS convergence history 
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The histogram for the bucklii g load P cr in Fig. 3 shows a normal probability distribution with a 
mean of 151,203 lb and a coefficient of variation of 2.8% (standard deviation of 4,245 lb). 
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Fig. 3 Histogram for axial buckling force, P„(lb) 

In examining the effect of the. coefficient of variation of the external force (i.e., CO V ( P a ) ) on 
cylinder failure probability, we found that in changing the COV from 1% to 10%, failure probability 
increases from Pj = 0.05 to Pf — 0.308, an increase of 516%. 

Although Monte Carlo simulation provides a convenient technique for the estimation of failure 
probability, it may be very inefficient when considering the computational time associated with each 
buckling analysis. A more efficient alternative to MCS involves the use of analytical techniques as 
discussed in next. 

Reliability Index 

If the limit state function is a linear function of all random variables, then we can use the first-order 
reliability method (FORM) to calculate the so called reliability index £ as the inverse of the 
coefficient of variation of the limit state function or simply 

tPcr-Vp (ID 

** I 9 2 

^P cr + <*P 
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where Up , Up are the mea 1 values of the critical axial buckling force and resultant axial force, 
respectively, and o P ,Up are the corresponding standard deviations. Since Pcr = Pcr(X),™ dall 
the random variables are assumed to be uncorrelated, we can obtain the variance of the buckling 
force using the partial derivative rule 



*k) 

dXi) 



( 12 ) 


If an analytical expression relating P cr to the pertinent random variables is available, then the partial 
derivatives in Eq. (12) can be obtained analytically, otherwise they may be found numerically using 
a finite difference scheme. 

If all random variables are normally distributed, then the probability of failure can be found as 
p f = i p(-p) ( 13 > 

where O is the cumulative distribution function of the standard normal variate. For example, 

P - 4.26 would correspond to Py = 10 5 . 


If, on the other hand, the linrit state function is a nonlinear function of random variables, then the 
method proposed by Hasofer and Lind 4 may be used to solve for fi . In this method, the original 
limit state, g(X) = 0, is transformed into the reduced limit state, g(X') = 0, with each reduced random 
variable found using the transformation 




Xj-Hx, 

a X, 


i = 1,2 ...,n 


(14) 


The point on the reduced limit-state surface g(X') = 0 having the maximum joint probability of 
failure is called the design point (also known as the most probable failure point, MPP) with 
coordinates identified by vector jc'*. The distance from the origin of the reduced coordinate system 

to the design point is the shortest distance to the failure surface, and is identified as /3or the 
Hasofer-Lind reliability index. Each coordinate of the design point is related to the reliability index 
as 



i 1, 2,. . ., 1 1 


(15) 


where a* represents the direction cosine corresponding to X ,, and is found as 



- , i,j - l,2,...,n 


(16) 


The asterisk symbol indicates that the derivatives of the reduced limit-state function are evaluated at 
the design point. Since the limit-state function is nonlinear, the solution for /3 becomes iterative as 
described in the appendix. 
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Reliability Sensitivity Analysi s 

A useful by product of analytical reliability technique is the evaluation of probabilistic sensitivity 
derivatives of (i with respect t .> the mean /I*, and standard deviation a Xj of each random variable 

as 


JP_ 

fax, 




i = 1,2 


(17) 


*2 

dp - ii 

d °x, P°x, (18) 

* ,* 

_ x i /—I? i 

°Xt 

Since the random variables (see Table 4) are of vastly different scales, the sensitivity derivatives are 
normalized as 




dP 


fyx, 


P 




dp 


da 


x, V 


P 


(19) 

( 20 ) 


To calculate /3 and its sensitivity derivatives, we must be able to evaluate the derivatives of the limit- 
state function with respect to individual random variables. These calculations, whether the limit- 
state function is linear or nonlinear, become very time consuming if there are many random 
variables and if g(X) is an implicit function of these random variables. 

To alleviate this burden and facilitate the subsequent reliability-based design optimization, we 
developed an analytical mod .4 of the limit-state function using the response surface methodology as 
discussed next. 

Response Surface Methodol ogy 

Response surface methodology (RSM) is mainly a statistical procedure used to develop a smooth 
mathematical function that represents an accurate functional relationship between the response 
variable and the independent parameters that influence it. 

A quadratic response surface model for n independent variables can be expressed as 


n n i 

AX, a) = ao + £ aft + £ £ aijXftj 


1=1 i=lj = 1 


( 21 ) 
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where oq, a t , a y are the unknown regression parameters for a total of (n + 1)(« + 2) / 2 including 
the interaction terms X t Xj, : * j. The accurate estimation of these regression parameters usually 
requires a large number of response samples, which could significantly increase the computational 
cost of the analysis. 

For estimation of cylinder reliability, we replaced the exact limit-state function in Eq. (7) by a close 
approximation in the form 

g(X)=P cr -P a (22) 

where P cr is the estimated ( itted) axial buckling force found using the algebraic response surface 
model defined by Eq. (2 1 ). 

For the cylinder specimen described in Table 4, we estimated the coefficients oq, a , fly by 
performing a least-squares l it of the regression model in Eq. (21) to the response values obtained 
from the Monte Carlo simulation. 

To validate the model, we examined various statistics such as the coefficient of determination R 2 and 
the root mean square error ( RMSE) found as 


RMSE = 



(23) 


where N is the number of response observations and M is the number of unknown coefficients in 
the response surface model. 

Since the approximate limit -state function in Eq. (22) is an explicit function of all 39 random 
variables, we would be able to calculate P and its probabilistic sensitivities through analytical 
differentiation of the response surface model. The mean buckling force estmated by the second __ 
order response surface equation is found to be 151,203 lb with model RMSE = 150.85 lb and R - 
0.9989. 


Table 5 below shows the values of the sensitivity derivatives of ft with respect to the mean and 
standard deviation of each random variable. The values in columns 3 and 5 are the normalized or 
logarithmic sensitivity derivatives that are normalized once more with respect to the largest 
sensitivity value in each group. 

As highlighted in Table 5, cf all random variables, the applied load. Young's modulus in the fiber 
direction, and cylinder diameter are found to have the greatest influence on cylinder reliability with 
axial buckling as the only failure mode. Although total shell thickness does have a large impac on 
cylinder buckling, the influence of each individual ply thickness is found to be less significant, it is 
also important to point out that the probabilistic sensitivity derivatives of ^ are greatly influenced by 
the variance in each random variable both implicitly and explicitly as indicated by Eqs. (17) and 
(18). 
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Table 5. Probabilistic sensitivity derivatives of reliability index (5 for the specimen in Table 4 


Random 

Variable 

Wlfrx ,) 

) norm 

(dp/ da Xi ) 

( ) norm 

E t 

8.14E-07 

8.44E-01 

-3.53E-07 

-2.90E-01 

e 2 

1 .68E-06 

1.30E-01 

-1.50E-07 

-1.24E-02 

V, 2 

6.92E-01 

9.01E-03 

-5.03E-03 

-8.1 IE-05 

0 I 2 

2.32E-06 

7.87E-02 

-1.65E-07 

-6.60E-03 

t, 

4.13E+02 

1.19E-01 

-7.93E+00 

-5.66E-04 

u 

4.40E+02 

1.27E-01 

-9.01E+00 

-6.43E-04 

2 

4.02E+02 

1.16E-01 

-7.49E+00 

-5.35E-04 

3 

L 

4. 1 0E+02 

1.18E-01 

-7.81E+00 

-5.58E-04 

4 

t. 

4.05E+02 

1.17E-01 

- 7.62E+00 

-5.44E-04 

5 

L 

3.93E+02 

1.13E-01 

-7.17E+00 

-5.12E-04 

6 

t 7 

4.42E+02 

1.27E-01 

-9.06E+00 

-6.47E-04 

to 

5.80E+02 

1.67E-01 

-1.56E+01 

-1.11E-03 

to 

5.63E+02 

1.62E-01 

-1.47E+01 

-1.05E-03 

t.ft 

4.. i i0E+02 

1.30E-01 

-9.40E+00 

-6.71E-04 

10 

t n 

3.46E+02 

9.95E-02 

-5.55E+00 

-3.96E-04 

1 1 

t,o 

3.38E+02 

9.73E-02 

-5.31E+00 

-3.79E-04 

I 2 

t.. 

3.42E+02 

9.85E-02 

-5.44E+00 

-3.88E-04 

1 3 

t.. 

3.52E+02 

1.03E-01 

-5.75E+00 

-4.10E-04 

1 4 

t,. 

5.28E+02 

1.52E-01 

-1.29E+01 

-9.23E-04 

t, 6 

5.12E+02 

1.47E-01 

-1.22E+01 

-8.70E-04 

0, 

-3 33E-02 

-8.63E-02 

-4.63E-04 

-2.98E-04 

02 

-2 79E-02 

7.22E-02 

-3.25E-04 

-2.09E-04 

0 3 

2.53E-02 

1.31E-01 

-5.36E-04 

-6.88E-04 

04 

1.86E-02 

9.61E-02 

-2.88E-04 

-3.70E-04 

e 5 

1.27E-02 

6.61E-02 

-1.36E-04 

-1.75E-04 

06 

7.92E-03 

4.10E-02 

-5.24E-05 

-6.74E-05 

07 

1.08E-02 

-2.79E-02 

-4.84E-05 

-3.1 IE-05 

08 

-4.49E-02 

-1.16E-01 

-8.43E-04 

-5.42E-04 

09 

-4.42E-02 

-1.15E-01 

-8.19E-04 

-5.25E-04 

0,o 

2.44E-02 

-6.32E-02 

-2.48E-04 

-1.60E-04 

e„ 

-2.43E-03 

-1.26E-02 

-4.95E-06 

-6.36E-06 

0,2 

-2.50E-03 

-1.29E-02 

-5.21E-06 

-6.69E-06 

0,3 

-1.43E-03 

-7.42E-03 

-1.71E-06 

-2.20E-06 

6,4 

-6.82E-05 

-3.54E-04 

-4.05E-09 

-5.20E-09 

0,5 

3 95E-02 

-1.02E-01 

-6.52E-04 

-4.19E-04 

0,6 

-4.35E-02 

-1.13E-01 

-7.90E-04 

-5.07E-04 

D 

2 85E-01 

2.46E-01 

-1.13E-02 

-2.43E-03 .... 

L 

-1.85E-01 

-1.49E-01 

-4.46E-03 

-8.92E-04 


-1.21E-04 

-1.00E+00 

-9.75E-05 

-1.00E+00 
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Figure 4 gives the plots of normalized sensitivity derivatives of (3 with respect to mean and standard 
deviation of each random vanable. 
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(b) 

Fig. 4 Normalized probabilistic sensitivities of ft with respect to (a) mean value and (b) standard 
deviation of each random variable for the specimen defined in Table 4 
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Effects of Distribution and efficient of Variation of Applied Load on J3 


Since the results of reliability analysis indicated a large sensitivity to the applied load, we also 
examined the effects of the distribution type and coefficient of variation of P a on (i and probability 
of failure. Two different distribution types (Normal and Lognormal) and three different values for 
coefficient of variation of P a were considered. In both cases, the applied load was assumed to have 
a mean of 143,690 lb. The critical buckling load was determined from the response surface model 
and was not affected by changes in COV of the applied load. 

The results shown in Table 6 below indicate that the coefficient of variation has a significant 
influence on both and probability of failure whereas the effect of the distribution type is relatively 
insignificant. It must be noted that in both cases the remaining 38 random variables were assumed 
to have normal distributions as specified in Table 4. 

T able 6. The effects o f d istribution and coefficient of variation of P„ on cylinder reliabili ty 

Reliability Index Method Monte Carlo Simulation 


P n Distribution COV (%) p P f 


Normal 

1 

1.7315 

0.042 

0.050 

5 

0.9099 

0.181 

0.186 


10 

0.5038 

0.309 

0.308 


1 

1.7326 

0.042 

0.050 

Lognormal 

5 

0.9188 

0.179 

0.184 

10 

0.5421 

0.295 

0.295 


Reliability-Based Structural Optimization 

The optimization problem is formulated as an element reliability index based structural optimization 
problem with a single design constraint imposed on buckling reliability index. Mathematical y, the 
optimization problem may be formulated as 

minimize W(X) 

subject to: > /3 min ^ 

Yf<Yi<Y?, i = 1,2,..., NDV 

where Wis the cylinder weight and f5 b is the reliability index associated with axial buckling. The 
design variables represent the mean values of only a subset of random variables such that Y, = (J-Xj 

with the corresponding standard deviations held constant during the optimization process. Each 
design variable is limited by lower and upper bound side constraints as indicated in Eq. (24). 

We chose a laminate concep. similar to that in Table 4, but with midplane symmetry (i.e., 

[±45 / 9 O 4 / +45] ), which reduces the total number of random variables from 39 to 23. Of the total 

of 23, only the 8 corresponding to ply thicknesses are treated as design variables. More 
specifically, their mean values are allowed to change while their standard deviations are heldtixed. 
The mean thickness for each ply is limited by lower and upper bounds of 0.0026 in. and 0.007 in., 
respectively. 
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The other geometric parame ers and material properties (in Table 4) are treated as random but with 
their means and coefficients of variation kept fixed during the optimization process. 

The optimization problem ir Eq. (24) is solved using the method of sequential quadratic 
programming in DOT 5 optimization software based on a global response surface modeling of the 
buckling response as explained next. 

Global Response Surface Te chnique 

In this case, a quadratic response surface model similar to that in Eq. (21) is developed for use over 
the entire design space. The population responses are generated using a direct Monte Carlo 
simulation with each random variable having a uniform distribution defined by the limits 

X],X l j=( l±g Xj )n Xj , j = 1,2,..., H (25) 

where g X j represents the bound increment for random variable Xj as defined in Table 7 below. 

By assuming each random variable is uniformly distributed, we are attempting to obtain samples 
from various points in the design space. This modeling of random variables, however, would not 
prevent the use of the resulting response surface model for design optimization where individual 
random variables are assumed to have a normal distribution. 

Table 7. Mean values and bound increments of random variables in Monte Carlo simulation 


Random Variable (Xp fi x . Sxj (%) 


Elast c Modulus, E t 
Elast c Modulus, E 2 
Poisson Ratio, v 12 

Shear Modulus, G 12 
Ply Thickness, r ply 
Ply Angle, 0 ply 
Cylinder Diameter, D 
Cyli n der Length, L 


18.0e6 psi 

4 

1.35e6 psi 

5 

0.226 

6 

0.543e6 psi 

6 

0.005 in. 

50 

45°, -45°, or 90° 

26 

15 in. 

50 

15 in. 

50 


The large increment on ply thickness allows for sampling over the entire feasible design space 
consistent with the corresponding lower and upper bounds specified in the optimization problem. 
We also chose large increments for cylinder diameter and length so the same response surface 
model can be used for optimization of cylinders with length and diameter of 10 to 20 in. Lastly, 
the support condition at the loaded edges is treated as deterministic with the external load assumed 
to be uniformly distributed over the top and bottom edges of the cylinder. 

We generated the necessary response data by conducting a total of 3,588 Monte Carlo simulation 
cycles. The buckling response for each cycle is found using the shell code with each interpolation 
function modeled by a 12 th degree Legendre polynomial. The buckling load found directly from 
this computer program is dubbed herein as "exact". 

The response surface model fitted through the sample data is found to have R 2 = 0 9774, COV = 
3.04% and RMSE = 4,383.26 lb. To assure the accuracy of the response surface model, we 
considered 20 random samples and compared the fitted buckling load to the exact value for each 
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sample. The results, shown in Table 8 below, indicate a maximum error of 4.22%, which is 
considered acceptable. 


Tabled . Comparison of exact and fitted buckling loads 

Random 

Sample 

Exact 

p„,ib 

Fitted 

P cr ,lb 

% Error 

1 

148551 

154305 

3.87 

2 

164791 

164622 

0.10 

3 

157294 

156206 

0.69 

4 

155613 

157132 

0.98 

5 

146943 

148753 

1.23 

6 

154598 

153307 

0.84 

7 

147618 

150250 

1.78 

8 

157701 

153539 

2.64 

9 

149712 

151621 

1.28 

10 

152033 

150656 

0.91 

11 

138397 

136995 

1.01 

12 

149330 

148239 

0.73 

13 

149616 

152113 

1.67 

14 

124890 

130162 

4.22 

15 

138740 

139190 

0.32 

16 

157662 

159231 

1.00 

17 

160469 

160284 

0.12 

18 

148128 

145619 

1.69 

19 

154835 

156182 

0.87 

20 

143155 

146534 

2.36 


Having developed the necessary response surface model, the optimization problem is s °' ve ^ 
the mean and coefficient of \ uriation of each random variable as specified in Table 9. I he solution 
to the optimization problem s given in Table 10 for three combinations of length and diameter and 
two different values for coefficient of variation of the applied load. The minimum reliability ind x 
is chosen as 4.26. The total shell thickness, h is shown in Table 10 in lieu of each mean ply 
thickness value. 

In all cases, the increase in COV (PJ causes an increase in the optimal laminate thickness an< J 
cylinder weight. The combi lation of small diameter and long length increases the buckling strengt 

of the cylinder. The large value for resulted in many of the thickness design variables to be 

pushed closer to the upper bound of 0.007 in. 


14 



Table 9. Statistical pr operties of random variables used in cylinder design optimization 


Fandom Dist. Mean COV 

Variable Type (%)_ 


p a 

Normal 

143,690 lb 

10 

E , 

Normal 

18.0e6 psi 

3.19 

E 2 

Normal 

1.35e6 psi 

4.26 

V 1 2 

Normal 

0.226 

5 

G 12 

Normal 

0.543e6 psi 

5.16 

? ply 

Normal 

Py, 

5 

0 Dlv 

Normal 

[±45/90 4 /+45] s 

5“ 

B 

Normal 

10, 15, 20 in. 

10 

L 

Normal 

10, 15, 20 in. 

10 


“Standard Deviation 


Table 10. Optimization results for (3 mm = 4.26 


'urameter 

COV = 5% 

COV = 10% 


fX D = 20 in., fi L = 

10 in. 

P cr , lb 

200,278 

223,995 

W, lb 

3,323 

3.542 

h, in 

0.0924 

0.0986 

Pd = Pl = 15 in 

P Cf ,lb 

194,019 

220,142 

W, lb 

3.692 

3.956 

h, in 

0.0912 

0.0976 


fi D = 10 in., fi L = 

20 in. 

Pen lb 

210,187 

230,097 

W,\b 

3.601 

3.790 

h, in 

0.0996 

0.1046 


Computational Requiremen ts 

Although the use of response surface model results in a considerable computational efficiency, it 
still requires further improvement. In this investigation, the major cost of the reliability-based 
design optimization was associated with the Monte Carlo simulation to generate the necessary poo 
of response samples for developing a reasonably accurate response surface model ot the buckling 
load. Once that part of the problem was completed, the actual optimization analysis took no more 
than one or two minutes. 

The 3,588 simulation cycles took approximately 448.5 CPU hours on a Sun server. Therefore, it is 
clear that other techniques for further improvement in computational efficiency must be 
investigated. 
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Appendix 

* * Jfc + 

For the composite cylinder, the nonlinear limit state function is expressed as g(.*i ,X2 ,---x n ) 

where (jc(*,X 2 *,.. x^) represent the reduced coordinates of the design point (or MPP), which are 
unknown initially, and arc calculated according to the iterative process outlined in Fig. A-l. 



Fig. A-l Piocedure for calculating the Hasofer-Lind reliability index 


For a more thorough descr ption of the procedure, the reader is referred to Refs. 4 and 6. 


17 







